void Navier::CompFM1(int i, int j, int k)//calculate temporary U velocity
{
                  int border=0;
                  double c,d,e0,e11,e12,f0,f11,f12,a1M1,b1M1;
                  double d2udx2, d2udy2, d2udz2;
                  double du2dxM1, duvdyM1, duwdzM1;
                  double Gx;
                  Gx = S.g[0]/S.froude;
                  //S.nu = 0.0;
                  a1M1=2*U[i][j][k] - U[i+1][j][k];
                  b1M1=U[i+1][j][k];
                  du2dxM1=DUU(a0,a1M1,U[i+1][j][k],b1M1,b0,DX[i],DX[i],DX[i],DX[i],S.alpha,border);
                    
                 if((j>=2) && (ISMOFLUID(flag[i][j-2][k])) && (ISMOFLUID(flag[i][j-1][k])))
                 {
                     e0=U[i+1][j-2][k];
                     e11=U[i+1][j-1][k];
                 }
                 else
                 {
                    if((j>=1) && (ISMOFLUID(flag[i][j-1][k])))
                    {
                        e11=U[i+1][j-1][k];
                    }
                    else
                    {
                        e11=U[i+1][j][k];
                    }
                  }
                  if((j<=S.gridp[1]-1) && (ISMOFLUID(flag[i][j+2][k])) && (ISMOFLUID(flag[i][j+1][k])))
                  {
                     f0=U[i+1][j+2][k];
                     f11=U[i+1][j+1][k];
                  }
                  else
                  {
                    if(j<=S.gridp[1]  && (ISMOFLUID(flag[i][j+1][k])))
                    {
                      f11=U[i+1][j+1][k];
                    }
                    else
                    {
                      f11=U[i+1][j][k];
                    }
                  }
                  if((i<=S.gridp[0]) && (j>=1) && (ISMOFLUID(flag[i+1][j-1][k])) && (ISMOFLUID(flag[i][j-1][k])))
                     c = V[i+1][j][k];
                  else
                     c = V[i][j][k];
                  if((i<=S.gridp[0]) && (j<=S.gridp[1]) && (ISMOFLUID(flag[i+1][j+1][k])) && (ISMOFLUID(flag[i][j+1][k])))
                     d = V[i+1][j+1][k];
                  else
                     d = V[i][j+1][k];
                  duvdyM1 = DUV(e0,e11,U[i+1][j][k],f11,f0,(DX[i+1]*V[i][j][k]+DX[i]*c)/(DX[i]+DX[i+1]),
                  (DX[i+1]*V[i][j+1][k]+DX[i]*d)/(DX[i]+DX[i+1]),DY[j-2],DY[j-1],DY[j],DY[j+1],DY[j+2],S.alpha,border);

                 if((k>=2) && (ISMOFLUID(flag[i][j][k-2])) && (ISMOFLUID(flag[i][j][k-1])))
                 {
                     e0=U[i+1][j][k-2];
                     e12=U[i+1][j][k-1];
                 }
                 else
                 {
                    if((k>=1) && (ISMOFLUID(flag[i][j][k-1])))
                    {
                       e12=U[i+1][j][k-1];
                    }
                    else
                    {
                       e12=U[i+1][j][k];
                    }
                 }
                 if((k<=S.gridp[2]-1) && (ISMOFLUID(flag[i][j][k+2])) && (ISMOFLUID(flag[i][j][k+1])))
                 {
                    f0=U[i+1][j][k+2];
                    f12=U[i+1][j][k+1];
                 }
                 else
                 {
                      if(k<=S.gridp[2]  && (ISMOFLUID(flag[i][j][k+1])))
                    {
                        f12=U[i+1][j][k+1];
                     }
                     else
                    {
                        f12=U[i+1][j][k];
                     }
                 }
                  if((i<=S.gridp[0]) && (k>=1) && (ISMOFLUID(flag[i+1][j][k-1])) && (ISMOFLUID(flag[i][j][k-1])))
                     c = W[i+1][j][k];
                  else
                     c = W[i][j][k];
                  if((i<=S.gridp[0]) && (k<=S.gridp[2]) && (ISMOFLUID(flag[i+1][j][k+1])) && (ISMOFLUID(flag[i][j][k+1])))
                     d = W[i+1][j][k+1];
                  else
                     d = W[i][j][k+1];
                  duwdzM1 = DUV(e0,e12,U[i][j][k],f12,f0,(DX[i+1]*W[i][j][k]+DX[i]*c)/(DX[i]+DX[i+1]),
                  (DX[i+1]*W[i][j][k+1]+DX[i]*d)/(DX[i]+DX[i+1]),DZ[k-2],DZ[k-1],DZ[k],DZ[k+1],DZ[k+2],S.alpha,border);

                  d2udx2=DD (U[i][j][k],U[i+1][j][k],b1M1,0,i);
                  d2udy2=DDS(e11,U[i+1][j][k],f11,1,j);
                  d2udz2=DDS(e12,U[i+1][j][k],f12,2,k);
                  F[i+1][j][k] = S.nu*(d2udx2+d2udy2+d2udz2) - du2dxM1-duvdyM1-duwdzM1 + Gx;
}